The onset of synchronization in large networks of coupled oscillators 
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We study the transition from incoherence to coherence in large networks of coupled phase oscil- 
lators. We present various approximations that describe the behavior of an appropriately defined 
order parameter past the transition, and generalize recent results for the critical coupling strength. 
We find that, under appropriate conditions, the coupling strength at which the transition occurs is 
determined by the largest eigenvalue of the adjacency matrix. We show how, with an additional 
assumption, a mean field approximation recently proposed is recovered from our results. We test our 
theory with numerical simulations, and find that it describes the transition when our assumptions 
are satisfied. We find that our theory describes the transition well in situations in which the mean 
field approximation fails. We study the finite size effects caused by nodes with small degree and 
find that they cause the critical coupling strength to increase. 
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I. INTRODUCTION 

In recent years, the importance of networks in different 
fields has become increasingly clear 0. [H Q ■ ^ has been 
observed that many real world networks possess topolo- 
gies which introduce important effects on the processes 
taking place on them. One of the most interesting and 
important of these processes is the synchronization of 
coupled dynamical systems. Synchronization is found in 
fields ranging from physics to biology 0, HJ , and in many 
cases involves a large network of dynamical systems. The 
structure of this network plays a crucial role in determin- 
ing the synchronization of the coupled elements. 

Kuramoto @ proposed and exactly solved a model for 
the synchronization of all-to-all uniformly coupled phase 
oscillators. His model and solution have become a guide 
as to how the coupling strength and the properties of the 
oscillators (e.g., their natural frequencies) might affect 
their synchronization, and generalizations of this basic 
model have been studied 0,13 . Some attempts to study 
the Kuramoto model with networks different from the 
all-to-all network have been made [8j . Networks in which 
the interaction strength depends on a distance have been 
studied, and it has been numerically found that a tran- 
sition from incoherent to coherent behavior occurs at a 
critical value of the coupling strength j^] ■ The Kuramoto 
model in networks without global coupling has recently 
started to receive attention. It was numerically observed 
[Tof that a transition is also present in scale free networks. 
Very recently, a mean field theory to determine the tran- 
sition to synchronization in more general networks has 



been proposed [Til IT^ | . The mean field theory result is 
that the critical coupling strength k m f is determined by 
the Kuramoto value, k , rescaled appropriately by the 
first two moments of the degree distribution of the nodes 
in the network: k m f — /cq (d}/(d 2 }, where 
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the degree d n of node n is the number of connections 
between node n and other nodes of the network, and N 
is the number of nodes in the network. 

In this paper we go beyond the mean field approxima- 
tion, obtaining a better estimate of the critical coupling 
strength. We also describe the behavior of a suitably de- 
fined order parameter past the transition. We show how 
our results reduce to those of the mean field theory when 
an additional assumption is introduced, and present ex- 
amples in different regimes. We find that in some regimes 
the mean field approximation does not provide an ad- 
equate description of the transition, whereas our more 
general estimate does. We also show how our results ex- 
plain observations for networks with distance dependent 
interaction strength. We study finite size effects caused 
mainly by nodes of small degree, and find that the tran- 
sition point is shifted to larger values of the coupling 
strength when these effects are taken into account. 

This paper is organized as follows. In Section [H] we 
present our theory and discuss the mean field approach. 
In Section IlIII we present numerical examples for differ- 
ent situations and test the different approximations. In 
Section llVI we discuss the case of networks with nonuni- 
form coupling strength. In Section^ we present a linear 
analysis of the problem. In Section IVII we consider finite 
size effects caused primarily by nodes with a small num- 
ber of connections. Finally, we conclude in Section IVIII 
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Some calculations were relegated to Appendices [XJ [El 
and O 

II. SELF CONSISTENT ANALYSIS 

As shown by Kuramoto [(|, the dynamics of weakly 
coupled, nearly identical limit cycle oscillators can, under 
certain conditions, be approximated by an equation for 
the phases O n of the form 

N 

On = W n + /J !Jnm(^m — On), (2) 
m— 1 

where uj n is the natural frequency of the oscillator n, N 
is the total number of oscillators and £l nm is a periodic 
function depending on the original equations of motion. 
The all-to-all Kuramoto model assumes that Q n m(O m — 
O n ) = (k/N) sm(O m — O n ), where k represents an overall 
coupling strength. In order to incorporate the presence 
of a heterogeneous network, we assume that 0, nm {O m — 
O n ) = kA nm sin(# m — O n ), where A nm are the elements of 
a NxN adjacency matrix A determining the connectivity 
of the network. Therefore, we study the system 

JV 

9„=w„ + i;^ A nm sin(6> m - O n ). (3) 

m— 1 

For specificity, we will primarily consider the case 
where the A nm are either (nodes n and m are not con- 
nected) or 1 (nodes n and m are connected, and all con- 
nections have equal strength). We assume that the net- 
work is undirected, so that A nm — A mn . We assume also 
that, for each n, the corresponding uj n is independently 
chosen from a known oscillation frequency probability 
distribution g(uj). We assume that g(uj) is symmetric 
about a single local maximum (cf . Sec [Vj , which with- 
out loss of generality we can take to be at w = 0. (If the 
mean frequency is loq ^ 0, we make the change of coordi- 
nates that shifts each w ra by ujq and each n by to^t.) In 
this case, synchronization will occur at frequency 0, i.e., 
n will remain approximately constant for synchronized 
nodes. 

We define a positive real valued local order parameter 
r n by 

r„e**» = A nm{e l9 ™)t, (4) 

m— 1 

where (•■•)< denotes a time average. In terms of r n , 
Eq. J2J can be rewritten as 

n = u n - kr n sm(O n - ij) n ) ~ kh n (t), (5) 

where the term h n (t) takes into account 
time fluctuations and is given by h n = 
Im{e- ie "J2m A nm((e ie ™) t - e* 8 ™)}, where Im stands 
for the imaginary part. Since we regard h n as a sum of 



d n approximately uncorrelated terms (where d n is the 
degree of node n given by d n — J2 m A nm), we expect h n 
to be of order \fd~n~. Substantially above the transition, 
due to the synchronization of the phases, the quantity 
r n w "}Z m A nm (e l9m )t is 0(d n ). Thus, if we assume that 
d n 3> 1, substantially above the transition the term h n 
can be neglected with respect to r n . However, just above 
the transition to coherence, the number of oscillators 
that are phase locked is small (see below), and so the 
term r n is also small. We need the number of locked 
oscillators to be large enough so that we can neglect 
h n , but, in cases where we use perturbative methods, 
we also require that the number of locked oscillators 
be small enough that the perturbative methods are 
still valid. We therefore do not expect the perturbative 
methods to agree perfectly just at the transition point. 
[Indeed in the classical Kuramoto (all-to-all) model a 
similar reservation holds for finite networks, as there are 

0(N^ 1 / 2 ) fluctuations of k YZ=1 gf "* for k below its 
critical transition value] In Sec. IVII we will investigate 
the effects of the time fluctuating term h n in Eq. 101, 
but, for now, we neglect it. 

With h n neglected in Eq. ®, oscillators with \uj n \ < 
kr n become locked, i.e., for these oscillators 9 n settles at 
a value for which 

sin(6>„ - ip n ) = u n /(kr n ). (6) 

(In general there are two such n ; the one closest to tp n 
is stable.) Then 

N 

r n =Y, A nm{e i ^-^ ) )t (7) 

m—1 



= V A e^ 9 --^) 

\u) m \>kr m 

In order to proceed further, we will introduce the follow- 
ing assumption: 

Assumption ~k We assume the existence of solutions 
Tn, ffin that are statistically independent of ui n . 

This is a nontrivial assumption; however, it is reasonable 
if most of node n's neighbors have reasonably large de- 
gree, so that they are not strongly affected by the value 
of bj n - And, as we show below, such a solution can be 
found in a self consistent manner. Using a milder version 
of Assumption -fc, we show in Appendix A that the sum 
over the unlocked oscillators in Eq. {7|) can be neglected. 
Therefore, only the locked oscillators remain in the sum, 
and we get from Eq. Q using Eq. ©, since r n is by 
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definition real, 



r n =Re{ J2 Anme^-^e^-^} (8) 

\uJm \ <krm 



= ^ A nm cos(tp m - i/) n )\ 1 



|w m |<fcr„ 



kr n 



V] A nm sin(?/> m - ip n ) ( ) , 



\u m \<kr 

where i?e represents the real part. For the imaginary 
part of Eq. (JZJ , we get 



= ^ A nm COsOm -1pn)[ 
\uj m \<~kr m 



V kr n 



(9) 



\<kr„ 



kTr. 



Using Assumption -fc, the contribution of the last term 
in the real part equation © can be neglected because 
of the symmetry of g(u) about 0. We thus obtain the 
approximation 



E 



An 



,COs(V>m-Vn)l/l- ( -j^— ) • (10) 



Since we are interested in the transition to coherence, 
we look for the solution of Eq. (|10|l that yields the small- 
est critical coupling k. The smallest critical coupling is 
obtained when the cosine in Eq. IjlOjl is 1. (Note that 
both the number of terms in the sum and their size de- 
creases as k decreases. Hence, a smaller k corresponds 
to a larger value of the cosine.) We therefore will look 
for solutions for which ip n — ip m = 0, i.e., ip n does not 
depend on n, and without loss of generality, we will take 
ijj n = 0. Note that this is a consistent condition in the 
sense that the imaginary part equation © is satisfied: 
the first term vanishes in the limit of a large number of 
connections per node due to the symmetry around of 
g(u>), and the second due to our assumed form that ipn 
does not depend on n. 

Equation (fTUfl then reduces to 



E 



A ,- 



kv r . 



(ii) 



If the particular collection of frequencies Lu n is known, 
this equation can be solved numerically. We will refer to 
this approximation, based on neglecting the time fluctu- 
ations in Eq. JSJ, as the time averaged theory (TAT). We 
now define an order parameter r by 



E 



iV 



(12) 



where al„ is the degree of node n defined by d. 
ELiAm- Note that r = E^=i d n {e^) t / ^ A 



coincides with the order parameter used in Refs. |lll 11^ 

If the number of connections per node is large, the 
particular collection of frequencies of the neighbors of a 
given node will likely be a faithful sample of the frequency 
distribution g(u>). Assuming this is the case, and using 
Assumption -^r, we approximate the sum in Eq. (|ll|l as 



r n — ^ y <A r 



— kr 



or, introducing z = uj/(kr m ), 

r n = ky^A nm r m / g{zkr m )\J\ - z 2 dz (14) 

This equation is one of our main results. It is analo- 
gous to Eq. (13) in Ref. [ill and E q- ( 6 ) in Ref - 
but, as opposed to including only information of the de- 
gree distribution of the network, it depends on the ad- 
jacency matrix, which completely describes the topology 
of the network. Equation l|14|) determines implicitly the 
order parameter r as a function of the network A nm , the 
frequency distribution g(u>), and the coupling constant 
k. We will refer to this approximation as the frequency 
distribution approximation (FDA). As with the TAT ap- 
proximation 1)11(1 . nonlinear matrix equation 1(14|1 can be 
solved numerically and the order parameter r computed 
from r n using Eq. i|12|) . 

We will now study the implications of Eq. (f 1 4|> by using 
approximation schemes in different regimes in order to 
obtain explicit expressions for the order parameter and 
the critical coupling strength. 

A. Perturbation Theory (PT) 

From the discussion above, coherent behavior is char- 
acterized by a nonzero value of r n . We determine the 
critical value of k by letting r„ — > + . The first order 
approximation g{zkr m ) f» <?(0) in Eq. (|14l) produces 



,(o) = ± 



/.•oE-' 1 '-" • 



,(0) 



(15) 



where fco = 2/(ng(0)). Since we are interested in the 
transition to coherence, the smallest k satisfying Eq. i|15|) 
is of interest. We thus identify the critical transition 
value of ko/k with the largest eigenvalue A of the adja- 
cency matrix A, obtaining 



kr — 



A 



(16) 



(In the case A nm = 1 of all-to-all coupling, A = N — 1.) 
Also rf is proportional to the mth component of the 



eigenvector u — [u±,U2, 



, u N 



T associated with this 
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eigenvalue. Note that this is consistent with Assump- 
tion -A-, since r n depends only on network properties (i.e., 
the matrix A) and is thus independent of oj n . Equa- 
tion (|16fl is one of our main results. It determines when 
the transition to coherence occurs in terms of the largest 
eigenvalue A of the adjacency matrix A. 

In order to assess how the order parameter r given by 
Eq. (|12[1 grows as k grows from k c , we must take into 
account that g(zkr m ) in Eq. i|14|) is not constant. For 
kr n small (see the discussion at the end of Sec. Ill B|) . the 
second order approximation yields 



Eqs. 112211 and (112.il) describe the behavior of the order 
parameter near the transition in terms of A and its asso- 
ciated eigenvector. We will refer to them as the pertur- 
bation theory (PT). 

The presence of the term (u 4 ) in Eq. 112MI) suggests 
that the expansion of g to second order might fail when 
there are a few components of the eigenvector u that are 
much larger than the rest. This occurs when the degree 
distribution is highly heterogeneous. We formulate more 
precisely this constraint in the discussion at the end of 
Sec. Ell 



(17) 

rn 

x J ^ (g(Q) + i 5 "(0)(zfcr m ) 2 ) v 7 ! - z 2 dz. 



Defining a = —irg"(0)ko/16, we get 
k 



k r X 



E A nm (r m - ak 2 

m } 



(18) 



We consider perturbations from the first order critical 
values as follows: 



r n = r { " ] + 5r n 



(19) 



where Sr n -C r^ <C 1 as k — > k c . Inserting this into 
Eq. (|TH|l . and canceling terms of order r„ , the leading 
order terms remaining are 



5r„ 



k r X 



k c X 

k-k. 
k r X 



r (o) 

nm 1 m ' 



In order for Eq. I|2U|) to have a solution for 8r n , it must 
satisfy a solubility condition. This condition can be ob- 
tained by multiplying by 7*n ^ , summing over 71, using 
Eq. I|15|) and the assumed symmetry A nm — A mn , to 
obtain 



T C 



,(0h4 



T ( 



.(Oh 



ak 3 



(21) 



In terms of u, the normalized eigenvector of A associated 
with the eigenvalue A, the square of the order parameter 
r can be expressed as 



for k/k c > 1, where 



(u) 2 X 2 
N(d) 2 (u A 



(22) 



(23) 



B. Mean field theory (MF) 

In this section we describe an approximation that 
works in some regimes and has the advantage of greater 
analytical tractability. In this section we also recover 
some of the results in Refs. Here we assume 

that r n is proportional to d n , r n cx d n . The assumption 
consists in treating the average 



N 



5> 

rn — 1 



(24) 



which depends on n, as if it were a constant independent 
of n. Following Refs. [TlTIT^. we call this the mean field 
(MF) approximation. It is also equivalent, near the tran- 
sition, to assuming that the eigenvector associated with 
the largest eigenvalue A satisfies u n cx d n . We will dis- 
cuss later the range of validity of this assumption. Note 
that this form for r„ is again consistent with our As- 
sumption if that r n is independent of uj n . The ratio 
r n /d n coincides under this approximation with the order 
parameter r defined in Eq. (|12|) . 

Summing over n and substituting r n = rd n in Eq. 114|l , 
we obtain 



N 



E d - 

rn — 1 



N 

k E 4 

m— 1 



g(zkrd m )\/ 1 — z 2 dz, (25) 



which coincides with Eq. (13) in Ref. ^lj. As we ap- 
proach the transition from above, r — > + , the first order 
approximation is g(zkrd m ) ~ <?(0), from which we obtain 



k — k m f — fcg 



\d 2 y 



(26) 



the main result of Ref. • 

In the limit N — > oo, we can replace (d q ) as defined by 
Eq. C) by 



(d q )oo = / d«p(d)dd, 



(27) 



where p(d) is the probability distribution function for 
the degree. Note that from Eq. (d q ) is always 

well-defined for finite N, but that Eq. I|27|) indicates 
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that (d q )oo diverges for power law degree distributions 
p(d) oc d~ 7 if 7 < q + 1. We also note that many real 
networks have approximate power law p{d) with 7 < 3 
(see Ref. p]). On the basis that (d 2 ) 00/(^)00 = 00 for 
2 < 7 < 3, Ichinomiya notes that from Eq. I|26() 
k m f — > as TV — > 00; i.e., predicts that in the limit 
TV — > 00 there is no threshold for coherent oscillations 
when 2 < 7 < 3. As will become evident, our numerical 
experiments, although for N 1, are often not well- 
approximated by the N — > 00 limit, in particular for 
7 < 3. 

The mean field approximation can be pushed fur- 
ther to second order by expanding g(zkrd m ) « g(0) + 
\g'' '(0)(zkrd m ) 2 in Eq. (|2~5)l . obtaining, provided Axc2 m is 
small, 



Z^m=l 



+ *V -«/'(()) 



(28) 



so that 



r - = [ ) ( _A_ _ 1 



for fc/fc c > 1, where 

(d 2 ) 3 



-^ 7— ' 7— (29) 



fe 



'/2 



(d*)(d) S 



(30) 



In expanding g to second order, it was assumed that kd m 
is small. The term (d 4 ) in Eq. (|30l) suggests that the 
conditions under which the expansion of g is appropri- 
ate are those under which (d 4 )^ is finite. In fact, Lee 
shows |l2j that for a power law distribution of the de- 
grees, p(d) tx d~ J , the above expansion is appropriate 

for 7 > 5. For 3 < 7 < 5, he obtains in the limit N — > 00 

/ , \ 1/(7-3) 
that r scales near the transition as r oc ( — 1 J 

A similar situation occurs in the perturbation theory 
[Eqs. (1221) and (|23[)]. which was also based on expanding 
g to second order. According to the previous discussion, 
we will only use the expression for r obtained from the 
perturbation theory for situations in which (d 4 )oo is fi- 
nite. The critical coupling strength in Eq. (|l(jf) . on the 
other hand, does not have this restriction. 

The expressions in Eqs. (|23|l and 13Uf) can be shown to 
coincide under the approximation u n oc d n . The treat- 
ment in Section [il Al does not assume that r n /d n is inde- 
pendent of n, and we will show in Section lllll that there 
are significant cases where it gives better results for the 
critical coupling strength than the mean field approxima- 
tion. 



C. Summary of approximations and range of 
validity 

In the previous sections, we developed different ap- 
proximations to find the critical coupling constant and 



Approximation 


Abbreviation 


Equation 


Time averaged theory 


TAT 


T^n 

11111 


Frequency distribution 
approximation 


FDA 




Perturbation theory 


PT 




Mean field theory 


MF 





TABLE I: Approximations considered, their abbreviation, 
and their corresponding equations. 



the behavior of the order parameter past the transition. 
Here we summarize the different approximations and the 
assumptions used in obtaining them. All the approxima- 
tions mentioned above assume that the number of con- 
nections per node is very large. This allowed us, among 
other things, to neglect the time fluctuating term h n (t) in 
Eq. (JSJ) ■ We will discuss the effect of this term in Section 

M 

The most fundamental approximation is given by 
Eq. Ullfl . This equation can be solved numerically if 
the frequency of each oscillator and the adjacency ma- 
trix is known. This is the time averaged theory (TAT). 
Assuming that the local mean field r n is statistically in- 
dependent of the frequency u n , the frequency distribu- 
tion approximation (FDA) given by Eq. (|14() is obtained. 
This equation can also be solved numerically, but only 
knowledge of the probability distribution for the frequen- 
cies and the adjacency matrix is required. Obtained by 
expanding the FDA approximation near the transition 
point, the perturbation theory (PT) describes the behav- 
ior of the order parameter in terms of the largest eigen- 
value of the adjacency matrix and its associated eigen- 
vector in networks where the degree distribution is rela- 
tively homogeneous, more precisely when (d 4 )oo is finite. 
Taking r n in the FDA approximation to be proportional 
to the degree, r n oc d n , leads to the mean field theory 
(MF). Table [I] summarizes the different approximations, 
their abbreviations and their corresponding equations. 
The diagram in Fig. ^ indicates the assumptions leading 
to each approximation. 

The mean field theory requires only knowledge of the 
frequency distribution and the degree distribution of the 
network, and thus it requires less information than the 
other approximations. However, it can produce mislead- 
ing results if not used carefully. The mean field approx- 
imation has the added assumption that the eigenvector 
u of A associated with the largest eigenvalue A satis- 
fies u n oc d n (since, close to the transition, r„ « u n ). 
While correlations might exist 13j, these two quanti- 
ties are in general not proportional. Further, the mean 
field approximation implies that A ~ (c? 2 )/(d), a re- 
sult that, although a good approximation in some cases, 
is not always true. Asymptotic forms for the largest 
eigenvalue in random networks with given degree dis- 
tributions are discussed and a sufficient condition for 
A ~ (d 2 )/(d) to be valid is presented in Q as follows. 
Let d m ax be the maximum expected degree of the net- 
work. If (d 2 )/(d) > y/d^\ogN, then A w (d 2 )/(d) 
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TAT Eq. (11) 




oUU 


co n distributed as g(co) 
independently of r n 






600 
400 




FDA 


Eq.(14) 




o n n 

zOO 
1 


k ~ k c , 
<d 4 > finite ^/ 








PT Eqs. (22,23) 






MF Eq. (25) 





< d >/<d> 




FIG. 1: Different approximations and the assumptions leading 
to them. See text for details. 
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FIG. 2: Largest eigenvalue A (diamonds) and (rf 2 )/(d) (stars) 
as a function of 7 for N = 5000 and do = 20. 



almost surely as N — > 00. We note also that, if the de- 
gree distribution is tightly distributed around its mean, 
so that \JT<P) ~ (d) ~ d max > (logiV) 2 , the condi- 
tion for the validity of A k (d 2 ) / (d) is satisfied. If in- 
stead \Jd max > ({d 2 } / (d))(log N) 2 , then almost surely 
the largest eigenvalue is A w yt as TV —> 00 01 • We 
will show that, indeed, to the extent that the approxi- 
mation A w (d 2 ) / (d) does not hold, the results from the 
numerical simulation of Eq. @ agree with the critical 
coupling strength as determined by the eigenvalue of the 
adjacency matrix, rather than by the quantity (d 2 }/(d). 

The asymptotic regimes described in 0] are not avail- 
able with the relatively small networks (N ~ 5000) we 
are restricted to study due to limited computational re- 
sources (see the end of Appendix |B}| . Also, finite but 
large networks are also interesting from an applied point 
of view. Thus, we numerically compare both approx- 
imations in order to illustrate the possible discrepan- 
cies between them in particular cases. Figure |21 was ob- 
tained using (for each 7) a single random realization of 
a network where the degrees d n are drawn from a power 
law degree distribution with power law exponent 7 (with 
d n > d = 20) and with N = 5000 nodes (see Sec. Iffl] 
for details on how the networks are generated). We plot 
(d 2 )/(d) and A as a function of 7. For the parameters 
used in the plot, (d 2 )/(d) coincides with the largest eigen- 
value A for values of 7 greater than 3. This suggests that 
the mean field result for the critical coupling strength 
k m f is valid for N = 5000 and 7 > 3. This is con- 
sistent with our numerical experiments in Sec. IIIII Wc 
show in Appendix [Bj however, that for 7 > 3 the mean 
field approximation (d 2 ) / (d) underestimates A for suffi- 
ciently large N (too large for us to simulate). In fact, as 
N — ► 00, A diverges while (d 2 ) / (d) remains finite in the 
case 7 > 3. Thus, the critical coupling constant obtained 
from our theory approaches zero as N — > 00, while the 
one obtained from the mean field theory remains con- 
stant. This suggests that the few nodes with high degree 
are able, for large enough N, to synchronize the network, 



and that these nodes are not taken into account by the 
mean field theory. 

For 7 < 3, we observe from Fig. [2] that A is less than 
(d 2 )/(d) when N = 5000. Thus, in this range, the mean 
field theory predicts a transition for a coupling constant 
that is smaller than that predicted by the perturbative 
approach. In the next section we will show, for a numeri- 
cal example in this regime, that the transition occurs for 
a larger coupling than that predicted by the mean field 
theory. 



III. EXAMPLES 

In order to test the results in Sec. |nj we choose a 
distribution for the natural frequencies given by g(u>) = 
(3/4)(l - uj 2 ) for -1 < uj < 1 and g(v) = otherwise. 
In order to generate the network, we specify a degree 
distribution and we use the "configuration" model (e.g., 
Sec. 4.2.1 of Ref. yij and references therein) to generate 
a random network realization with the specified degree 
distribution: (i) we first generate a degree sequence by as- 
signing a degree d n to each node n according to the given 
distribution; (ii) imagining that each node n is given d n 
spokes sticking out of it, we choose pairs of spoke ends 
at random, and connect them. 

We consider a fixed number of nodes, N = 2000, and 
the following networks with uniform coupling strength 
(i.e., A nm — 1 or 0) (i) the degrees are uniformly dis- 
tributed between 50 and 149, and (ii) the probability of 
having a degree d is given by p(d) oc g?~ 7 if 50 < d < 2000 
and p(d) = otherwise, where 7 is taken to be 2, 2.5, 3 
and 4. [Our choice p(d) = for d < 50 insures that 
there are no nodes of small degree, and suggests that 
our approximation of neglecting the noise-like, fluctuat- 
ing quantity h n in Eq. J5J is valid. We return to this 
issue m 

Section EH] 

The initial conditions for Eq. J3J are chosen randomly 
in the interval [0, 27r] and Eq. © is integrated forward 
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in time until a stationary state is reached (stationary 
state here means stationary in a statistical sense, i.e. the 
solution might be time dependent but its statistical prop- 
erties remain constant in time) . From the values of 9 n (t) 
obtained for a given k, the order parameter r is esti- 



mated as 



Em=l^ m ( e ' (m ) f /Em=l 

time average is taken after the system reaches the sta- 
tionary state. (Close to the transition, the time needed 
to reach the stationary state is very long, so that it is dif- 
ficult to estimate the real value of r. This problem also 
exists in the classical Kuramoto all-to-all model.) The 
value of k is then increased and the system is allowed to 
relax to a stationary state, and the process is repeated 
for increasing values of k. 



where the 



1.0 
0.8 
0.6 
0.4 
0.2 



▲ 


Simulation 




TAT ^ 




PT 








MF /{ 











0.5 



1.5 



k/k, 



The results for the networks with power law degree dis- 
tributions [networks (ii)] are shown in Figs. 01(a), (b), (c) 
and (d) for 7 = 2, 2.5, 3, and 4, respectively. The order 
parameter r 2 from numerical solution of the full system 
in Eq. (J3J) (triangles), the time averaged theory (solid 
line), the frequency distribution approximation (stars), 
and the mean field theory (long-dashed line) are plotted 
as a function of k/k c . We do not show the perturbation 
theory since in all these cases 7 < 5 and so we do not 
expect the perturbative theory to be valid as N — » 00. 

The time averaged theory agrees best with the numer- 
ical simulations in all cases. The frequency distribution 
approximation also agrees well in all cases; though it pre- 
dicts a sharper transition than actually occurs. The mean 
field approximation agrees closely with the frequency dis- 
tribution approximation for 7 = 4 and, away from the 
transition, for 7 = 3. However, for 7 = 2 and 7 = 2.5, it 
deviates greatly from the other approximations and from 
the numerical simulation. The critical coupling strengths 
predicted by the mean field theory and by the perturba- 
tion theory are very close for 7 = 4, but the mean field 
theory predicts a transition at about 10% smaller cou- 
pling for 7 = 3, about 20% smaller for 7 = 2.5, and 
about 40% smaller for 7 = 2. Since the transition in the 
numerical simulation is not so well-defined, both approx- 
imations are reasonable for 7 = 3, but for 7 = 2 and 
7 = 2.5 the critical coupling strength predicted by the 
mean field approximation is clearly too small. 



FIG. 3: Order parameter r 2 obtained from numerical solu- 
tion of Eq. © (triangles), time averaged theory (solid line), 
mean field theory (long-dashed line), and perturbation theory 
(short-dashed line) as a function of k/k c for network (i), with 
the degree of the nodes uniformly distributed in {50, . . . , 149}. 
All curves are obtained using the same single random network 
realization. 



In Fig. |3 we show the results for the network with a 
uniform degree distribution as described above [network 
(i)]. We plot r 2 from numerical solution the full system 
in Eq. © (triangles), the theoretical prediction from the 
time averaged theory (solid line) , the prediction from the 
mean field theory (long-dashed line), and from the per- 
turbation theory (short-dashed line) (see Table as a 
function of k/k c , where k c is given by Eq. (|16() . The fre- 
quency distribution approximation agrees with the time 
averaged theory, so we do not include it in the plot. In 
this case, all the theoretical predictions provide good ap- 
proximations to the observed numerical results. The time 
averaged theory reproduces remarkably well the numeri- 
cal observations. Even the irregular behavior near the 
transition is taken into account by the time averaged 
theory. The mean field theory is in this case a good 
approximation, providing a fair description of the order 
parameter past the transition. The perturbation theory 
is valid in this case up to k/k c « 1.3. 



In the past years, it has been discovered that many real 
world networks have degree distributions which are power 
laws with exponents between 2 and 3.5 0,0 El- 111 01 '- 
der to accurately predict the critical coupling strength 
across this range of exponents, the critical coupling con- 
stant given by k c = ko/X determined by the largest eigen- 
value of the adjacency matrix should be used. The be- 
havior of the order parameter can be estimated using the 
time averaged theory or the frequency distribution ap- 
proximation. These two approximations were found to 
be consistently accurate for the range of exponents and 
values of the coupling constant studied. For the value 
of N used, the mean field theory works well in predict- 
ing the critical coupling strength and the behavior of the 
order parameter if one is interested in values of 7 larger 
than 3. 

Tables [n] and II I II present the results of comparing the 
theoretical predictions with the numerical integration of 
Eq. f° r different networks. Table [H] compares the 
observed critical coupling strength with the theoretical 
estimate. If both are close, the entry is "G" , and other- 
wise "NG". Table UTTl compares the predicted behavior of 
the order parameter past the transition with the observed 
one. If the corresponding entry in Table [H] is "NG", no 
comparison is attempted. The entries are the range of 
k/k c over which the corresponding theoretical prediction 
agrees with the numerical simulation. 
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FIG. 4: Order parameter r obtained from numerical solu- 
tion of Eq. (triangles), time averaged theory (solid line), 
frequency distribution approximation (stars), and mean field 
theory (long-dashed line) as a function of k/k c for degree dis- 
tributions given by p(d) oc d -7 if 50 < d < 2000 and p(d) = 
otherwise, with (a) 7 = 2, (b) 7 = 2.5, (c) 7 = 3, and (d) 
7 = 4. All curves in each figure are obtained using the same 
single random network realization. 



Degree distribution 


TAT 


r DA 


1\ TP 

Mr 


OT 

r 1 


p(d) uniform in {50, . . . , 149} 


G 


G 


G 


G 


p(d) oc cT 7 , 7 = 2 


G 


G 


NG 


G 


p(d) oc d~ 7 , 7 = 3 


G 


G 


G 


G 


p(d) oc d -7 , 7 = 4 


G 


G 


G 


G 



TABLE II: Comparison of the predicted critical coupling 
strength versus the observed one for the different approxima- 
tions (columns) and different networks (rows). If the critical 
coupling strength is predicted by a given approximation for 
a certain network, the corresponding entry is marked "G". 
Otherwise, "NG" is entered. 
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TABLE III: Comparison of the predicted behavior of the or- 
der parameter versus the observed one for the different ap- 
proximations (columns) and different networks (rows) . If the 
behavior is correctly predicted by a given approximation for a 
certain network, the corresponding entry contains the range of 
k/k c after k/k c = 1 for which the approximation works well. 
A "+" indicates that the agreement possibly persists for larger 
values of k. When "NG" appears in the corresponding entry 
in table [H] no comparison is attempted and a "-" is entered. 
A "-" is entered when the perturbation theory is inapplicable 
(7 < 5), see Sec. EH 



d n = ^2 m A nm , our results carry through to the more 
general case of non uniform coupling. As an example of 
this situation, we apply our results to the case treated 
in Refs. of a distance dependent interaction strength. 
Assume that the nodes n are equidistantly located on a 
circle and the matrix elements are given by 



A nm = f(\n - m|), 



(31) 



where \n — m\ represents distance modulo N (e.g. 
|1 - N\ = 1), /(0) = 0, and / > 0. Then each row of A 
has the same sum A = ^ m A nm , and [1, 1, ... , 1] T is an 
eigenvector with eigenvalue A. By the Gershgorin circle 
theorem (each eigenvalue a of A satisfies, for some 
n, \a- Ami < E ra ^„ lAiml), tms is tu e largest eigen- 
value (since A nn = 0), and thus determines the transition 
to synchrony as described in the previous section. This 
scaling factor has been proposed before, by analogy to 
spin systems, to determine the transition to coherence 
in the case of a power law decaying interaction strength 
/(x) = x" 7 0. 



IV. NONUNIFORM COUPLING STRENGTH 

So far, our examples have assumed that the coupling 
strength is uniform (i.e., all the entries of the adja- 
cency matrix A have been taken to be or 1). How- 
ever, considering that the degree of a node is defined as 



V. LINEAR STABILITY APPROACH 

Partly as a precursor to the next section (Sec. IVI|I . 
in this section we discuss another approach that has the 
advantage of providing information on the dynamics of 
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the system. We study the linear stability of the incoher- 
ent state by a method similar to that used in Ref. [l7j . 
We assume that in the incoherent state the solution to 
Eq. is given approximately by 



= LO n t + i 



(32) 



where <p n is a random initial condition. We introduce 
infinitesimal perturbations to this state by 



9« — n + $n 



(33) 



In Appendix El we assume that the perturbations grow 
as a function of time as e st , and obtain the eigenvalue 
equation 



k 



N 



A 



m— 1 



nm^m 



%UJ r , 



(34) 



We look for solutions b n of this equation that are inde- 
pendent of the frequencies io n (similar to Assumption -jf). 



Under this assumption, replacing (s — iu n ) 
with its expected value, we get 



K = \ \ s 



1 



N 

E 

m— 1 



Eq. 



(35) 



where 



g(u>)du> 



iuj 



(36) 



and the integration contour is defined in the causal sense 
[i.e., for Re(s) > it is along the real axis, and for 
Re(s) < it passes above the pole oj — —is]. We thus 
obtain the dispersion relation 



^ k\ f g(u>)duj 



(37) 



where, as in Sec. [HI A is the largest eigenvalue of the adja- 
cency matrix A. Except for the presence of the eigenvalue 
A, this is the known dispersion relation for the stability of 
the incoherent state of the Kuramoto model Q • Under 
our assumption that g(u>) is even and decreases mono- 
tonically away from (Sec. UJ, an unstable [Re(s) > 0] 
solution of Eq. H37|) is real (note that, since A is sym- 
metric, A is real). In order to find the critical coupling, 
we let s — ► + , (s — iuj)^ 1 — > iP(l/uS) + n8(u>). Since 
g(oj) is symmetric, ((s — iuj)^ 1 ^ — > ng(0). According to 
Eq. (J33J), the critical coupling is then given by 



VI. EFFECT OF FLUCTUATIONS 

So far we have neglected the effect of the small fluctu- 
ations due to the finite number of connections per node. 
In our examples, we have presented networks that do 
not have nodes with small degree. However, in many 
networks there is a large fraction of the nodes with small 
degree; in all our examples in Sec. lllll there were no nodes 
with degree less than 50 [p(d) = for d < 50]. For ex- 
ample, scale free networks generated using the Barabasi- 
Albert method |2| sometimes have parameters so that 
(d) = 6. 

In developing our theory, we neglected the time varia- 
tions in Eq. and worked thereafter with the average 
value of the phase of the locked oscillators. In order to 
gain insight into the effect of these fluctuations, we will 
treat the time fluctuations as a noise term. 

The theory we present is heuristic and may be thought 
of as an expansion giving a small lowest order correction 
to the linear stability approach of Sec. 0for large but 
finite (d). On the other hand, later in this section, we will 
apply this theory to numerical examples where the finite 
size effect is not small, and we will find that the theory is 
still useful in that it correctly indicates the trend of the 
numerical observations. 

Like in Sec. we consider perturbations to the in- 
coherent state described by Eq. I|32|l . As an approxi- 
mation, we regard the coupling term in Eq. f n (t) = 
fc^ =1 A nm s'm(9 m — 8 n ), as a noise term. In addition to 
growing linearly with time, the phase of the oscillator n 
will diffuse under the influence of this noise. We assume 
that 9 n {t) = (f> n +u) n t + W n (t), where W n (t) is a random 
walk such that (W n (t)) = and (W m (t)W n (t)) = 2D nm t, 
and 4> n is an initial condition, which we assume to be 
randomly drawn from [0,2tt). (In this section, by (. ..) 
we mean an expected value, i.e. and ensemble average, 
rather than an average over t or n.) 

By using the linear approach of Section the dif- 
fusion coefficients D nm will give us information on how 
the critical coupling strength differs from Eq. J3SJ. The 
diffusion coefficients D nm are given by 



D nm = / {f n (t + r/2)f m (t-r/2))dT (39) 



A nj (sm(8+ - 6+)A mk sin(0 fc - 6> m ))dr, 



k 

A ' 



(38) 



in agreement with the nonlinear approach. [Wc note, 
however, that, if g(u>) has multiple maxima, then the first 
instability can occur at Im(s) ^ at a value of k below 
that given by Eq. (|38p. This is why we have assumed 
that g(io) decreases monotonically away from u> = 0.] 



where + (respectively — ) indicates evaluation at t + r/2 
(respectively t—r/2). Consider first the case n ^ m. The 
contribution of the terms with {j, n} ^ {k, m} vanishes 
after the integration, and we obtain, using the symmetry 
of A, 



D„ 



9+) S m(6--6-)). (40) 
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We now introduce our aforementioned assumption that 
n (t) — u) n t is a random walk plus a random initial con- 
dition, 9 n (t) — <fi n + w n t + W n (t). Using the identity 
sin(a;) sin(y) = (cos(x — y) — cos(x + y))/2 and averaging 
over the initial phases <fi n we get 



Al m (cos(AW m - AW n + u; mn T))dT, (41) 



where AW n = W 7 f — W~ and ui mn = iv m — w n , We now 
use the fact that for a Gaussian random variable x with 
variance a 2 we have (cos(x)) = Re(e lx ) — Re(e i ( x '*~ a / 2 ). 
In our case, (x) = w mn r and a 2 = ((AW m — AW n ) 2 ) = 
2(D n + D m - 2D nm )T, where D n = D nn . After using 
this to compute the expected value, and performing the 
integration, we obtain for n =/= in 



Dn 



k 2 

— A 2 

2 m 



D n + D r . 



2D,, 



{Dr, 



D — 2D ^ 2 



(42) 



If n = to, the calculation proceeds along the same lines, 
but the nonvanishing terms in Eq. (|39|l are those for 
which k = j. Together with Eq. JUJ, this results in 



D„ 



E 



D„ 



(43) 



In principle, Eqs. Q42J) and (|43|) can be solved for D n 
as a function of k if the frequencies and the adjacency 
matrix are known. 

In order to relate the diffusion coefficients to the criti- 
cal coupling constant, we resort to the linear analysis of 
Sec.[V] When noise is introduced in the linear approach, 
Equation (|34|l for the growth rate s generalizes, as shown 
at the end of Appendix |0 t° 



i N 

k X ^ 

9. 2.*, 



Dr, 



lU> r , 



(44) 



Since Re(s) > corresponds to instability of the inco- 
herent state, it is expected that the effect of the noise as 
reflected by positive D m is to shift the transition point 
so that the critical coupling constant is larger. 

In order to solve for the growth rate s for a given value 
of k, we rewrite Eq. (|44|l as 



b = -D(s)Ab, 



(45) 



where b is the vector with components {&„}, D(s) is the 
diagonal matrix D(s) = diag{(s + D m — iw m ) _1 }, and A 
is the adjacency matrix. The characteristic equation is 



drt | '-D(s)A- 



I 



0. 



(46) 



where / is the N x N identity matrix. This implies 



<!<•: ( ^A-Dis)- 1 ) =(). 



(47) 



det ( —A — diag{D r , 



iujm} - si 



0. 



(48) 



that is, the growth rate s is an eigenvalue of the matrix 
M(k) = (k/2)A - diag{D rn - iuj m }- 

For a given value of k, Eqs. (|42|l and (|43|l can be solved 
iteratively. We have found that, by starting from an ini- 
tial guess for the values of D nm and repeatedly evaluat- 
ing the right hand side of Eq. (|42J) in order to get the 
next approximation to the values of D nm , convergence is 
achieved to a solution that is independent of the initial 
guess if the condition D n > is imposed. When the val- 
ues of D nrn have been found for a given value of k, the 
relevant growth rate is calculated as the largest real part 
of the eigenvalues of the matrix M(k) defined above. 

As an example, we consider three networks with the 
degree of all nodes d given by d = 100 in the first, d = 50 
in the second and d — 20 in the third one. In order to 
solve numerically the coupled equations, we work with a 
small number of nodes, N — 500. In Fig. [S] we show the 




FIG. 5: Order parameter r 2 obtained from numerical solution 
of Eq. (solid lines) and growth rate Re(s) (dashed lines) 
for a network with the degree of all nodes d = 20, d — 50 
and d = 100 as a function of k/k c . The arrows indicate which 
network corresponds to the given curve. 



results for a realization of the three networks. The order 
parameter r 2 obtained from numerical solution of Eq. @ 
(solid lines) and the growth rate obtained from Eqs. 142fl . 
(|^5Tl and (dashed lines) are plotted as a function of 
k/k c . The arrows indicate which network corresponds to 
the given curve. We observe that, as the connections per 
node are decreased, the transition point shifts to larger 
values of the coupling constant. This trend is reproduced 
by the growth rate curves, which are displaced to the 
right for smaller values of the degree. 

We emphasize that the theory we described above is 
applicable to networks for which (d) is large but finite. 
However, in Fig.[S]we applied the theory to cases in which 
(d) is not very large. Although we do not expect the 
theory to be valid in this case, we find that it correctly 
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describes the trend present in the numerical observations, 
i.e., a shifting of the transition to coherence to larger 
values of the critical coupling as nodes of small degree 
become important. 



VII. DISCUSSION 

A transition to coherence in large networks of cou- 
pled oscillators should be expected at a critical value of 
the coupling strength which is determined by the largest 
eigenvalue of the adjacency matrix of the network and its 
associated eigenvector. In the all-to-all case, the largest 
eigenvalue is N — 1 w N and thus the Kuramoto result 
k c = kci/N is recovered. The largest eigenvalue of the 
adjacency matrix of a network is of both theoretical and 
practical importance, and thus its properties have been 
studied in some detail [la.H^ . We remark that our anal- 
ysis allows the case of nonuniform interaction strengths 
by introducing continuous values in the entries of the ad- 
jacency matrix A. 

We developed different approximations in order to de- 
scribe the transition to coherence in terms of an appro- 
priately defined order parameter which generalizes the 
parameter used in the classical Kuramoto model |ll| . 
See Table H] and Figure ^ for a summary of the approx- 
imations and assumptions. The time averaged theory 
(TAT) provided the most accurate description of the be- 
havior of the order parameter, and assumes knowledge of 
the adjacency matrix A nm and the individual frequencies 
u) n . The frequency distribution approximation (FDA) 
also provides a good approximation but does not require 
knowledge of the individual frequencies. These approx- 
imations yield equations that have to be solved numeri- 
cally. The time required to numerically solve these equa- 
tions is, however, much less than that required to numer- 
ically integrate the original differential equations. The 
perturbation theory (PT) yields analytic expressions for 
the order parameter when close to the transition in terms 
of the largest eigenvalue of the adjacency matrix and its 
associated eigenvector, but is limited to networks with 
a relatively homo gene ous degree distribution. The mean 
field theory (MF) [i"H is obtained by introducing the ad- 
ditional assumption that the components of the eigen- 
vector associated with the largest eigenvalue are propor- 
tional to the degree of the corresponding node. This 
does not necessarily have to be the case when close to 
the transition, and because of this extra assumption, we 
expect the other approximations to more generally accu- 
rately describe the transition than the mean field theory. 
Figs. 0Ja) and (b) show that for the particular case of 
scale free networks with N = 2000, 7 = 2 and 7 = 2.5 
this is the case. In general, we observed that for low val- 
ues of the exponent 7 (see Fig.0} the mean field approx- 
imation and the perturbative approximation yield differ- 
ent critical coupling strengths. The mean field theory 
has the advantage that analytic expressions can be com- 
puted without the need of solving the eigenvalue problem 



for the adjacency matrix, and could be useful when only 
limited information is available about the network. How- 
ever, in general, our results suggest that one of the other 
approximations mentioned above should be used. 

We remark that even though the time averaged theory, 
the frequency distribution approximation and the pertur- 
bation theory require in principle knowledge of the full 
matrix A, knowledge of the degree distribution may be 
enough in some cases. As in our examples, an adjacency 
matrix A can be generated randomly with a given degree 
distribution. Our results indicate that even this limited 
reconstruction of the original network might improve the 
mean field results (see Sec. lIIII) . 

Our assumptions restrict the class of networks for 
which the results apply. We assumed that sufficiently 
near the onset of synchronization each node is coupled 
to many locked oscillators. In practice this implies that 
most nodes should have a high degree. This is an im- 
portant restriction for our theory. In Sec. II I II we used 
networks with a minimum degree of 50. As mentioned 
before, we observed that in networks with small aver- 
age degree (about 20) , the observed critical coupling was 
larger than the one predicted by our theory. By includ- 
ing the previously neglected time fluctuations, we devel- 
oped a heuristic theory in Sec. IV1I which correctly pre- 
dicts the trend observed in the numerical simulations. 
As the nodes with small degree become important, both 
our theory and the numerical observations indicate that 
the transition to synchrony occurs at larger values of the 
coupling strength. 

In conclusion, we have developed a theory predicting 
the critical coupling for the transition from incoherence 
to coherence in large networks of coupled oscillators. We 
found that for a large class of networks, a transition to co- 
herence should be expected at a critical value of the cou- 
pling strength which is determined by the largest eigen- 
value of the adjacency matrix of the network. We devel- 
oped and compared various approximations to the order 
parameter past the transition, and studied the effect of 
the fluctuations caused by finite size effects. 

This work was sponsored by ONR (Physics) and by 
NSF (contracts PHYS 0098632 and DMS 0104087). 

APPENDIX A 

In this Appendix we show that, using Assumption -fc, 
we can neglect the sum over the unlocked oscillators in 
Eq. 0, 

JV 

J2 A nm (e ie ™) t . (Al) 

\LJm I >kr rn 

We will follow to some extent Chapter 12 of Ref. [4]. The 
time average is given by 

(e^) t = r e ie p m (6)d0. (A2) 
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where p m (9)d9 is, given the connections of node m and its 
natural frequency u> m , the probability that its phase 9 m 

lies in the interval [9,6 + d9). It satisfies p m (9) oc 1/ 9 . 
Including the normalization we have, neglecting the term 
h m in Eq. (0, 



The sum in Eq. (|A1|) can be written as 



Pm{9) = 



k 2 r 2 



2tt \ui r . 



kr m sin(9 - ip m )\ ' 



(A3) 



a i ie m \ a /~~2 r~2~ ■ / \* r e t6 {uj m + kr m sm(9-ip m ))d9 



The integral of the first term vanishes since the 2ir- 9 — > 6* + 7r. We are left with 
periodic integrand changes sign under the transformation 



E X a i ie \ a r^r> / \ 1 C el8 sin(0 — ^ m )d9 . 
A nm (e m )t = > ^mV^-^mWrnh- -5 ;2 2 . 2 la r~ T ( A5 ) 



In this sum, sign(uj m ) is independent of li^ and, using 
Assumption -fc, it is independent of r„ and as well. If 
there are many terms in the sum, it will be then of order 
\/dn due to the symmetry of the frequency distribution, 
and thus will be small compared with the sum over the 
locked oscillators, which is of order d n [see Eq. <|11|) ]. 
Note that we did not use here the full strength of As- 
sumption since we only required the sign of u> m to be 
independent of r m and tp m . 



{d 2 } I (d) ~ do, where do is the minimum degree \p(d) = 
for d < do]. The maximum degree is approximately given 
by d max ~ doN 1 ^ 1-1 ^ [J. Inserting these estimates in 
the condition y 'd max ~ (d 2 ) / (d) (log N) 2 we obtain 

iV^-^logAT) 4 ^- 1 ). ( B1 ) 

As an example, for 7 = 4 and do = 20, the upper bound 
is approximately N ~ 10 25 , a far larger system than we 
can simulate. 



APPENDIX B 

Here we show that for sufficiently large N and a power 
law degree distribution of the degrees, p(d) oc d~"> , the 
mean field approximation (d 2 ) / (d) underestimates A for 
7 > 3. We base our argument in the results of Ref. [l4| : 
if for a random graph y 'd max > (d 2 ) / (d) (log N) 2 , then 
max almost surely as — > oo, where d max is the 
largest expected degree. 

In the case under consideration (7 > 3), (d 2 )/(d) con- 
verges to the finite value (d 2 ) 00 /(d) 00 [(. . . ) O0 is defined 
by Eq. while d max diverges as JVVfr-i) |^. T hus, 

for large enough N, the conditions for A ~ sj d max will 
be satisfied, since — * 00 as N — > 00. 
While A ~ V d max —* 00 as N — + 00, the mean field 
approximation (d 2 )/(d) remains finite. 

We can estimate an upper bound on how large N needs 
to be for this discrepancy to be observed. For large N, 



APPENDIX C 

In this Appendix we study the linear stability of the 
incoherent state by a method similar to that presented 
in Ref. [T3. As described in Section [Vl we assume that 
in the incoherent state the solution to Eq. © is given 
approximately by 

0£(t)«Wn* + ^, (CI) 

where <fi n is an initial condition. We introduce infinitesi- 
mal perturbations to this state by 

9 n = 0° n + 6 n . (C2) 

Linearizing Eq. 1(3}, we get 

N 

5 n = A nm cos(9 m - 9°)6 m + //„ - v n 5 n , (C3) 

m— 1 
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where (i n = k E TO =i A nm sin(0^ - 0°) and v n = 
k ELi 4m cos (^m — ^n)- As before, we assume that 
the number of links to node n is so large that, due to the 



incoherence, we may neglect the terms \x n and v n . With 
this simplification, Eq. l|C3f) can be recast as an integral 
equation as follows: 



S n (t) = k ( dt'Y, A nm 8 m {t') cos[e Q m {t') - 9° n (t')} (C4) 

m=l 

1 ft ( N N \ 

"'- 00 \m=l m=l / 



Multiplying by A^e* 9 *^*), summing over n and defining 
Bn(t) = E™=1 A im <W*)e^ (t) , we get 



Si(t) = 



iv 

n=l 



(C5) 



We assume that the quantities B„ grow exponentially 
with time as B n (t) = b n e st , where Re(s) > 0. Inserting 
this ansatz in Eq. I|C5|) . and performing the integration 
we get 

h N A ■ h h N A- h*p 2i8 °W 

_ « \ - A Jn O n K 2ilm(s)t \ " ^""n 6 



&i = t: — " — +- -' 

n—1 



y^uni . (C6) 

n—1 



The second sum is very small due to the incoherence 
of the 0°'s. So, changing indices, we are left with the 
eigenvalue equation 



N 



k \ A. nrn b 

6 « = ;L 



nm w m 



2 z — ' 5 — 2CcJ r , 
m— 1 



(C7) 



as claimed in Section Ivl 

If, as proposed in Section IVT1 there are fluctuations in 
the values of d°(t) such that 9°(t) = uj n t + <j) n + W n (t), 
where W n (t) is a random walk such that (W n (t)) — and 
(W n (t) 2 ) = 2D n t, we take the expected value of Eq. IIC5I) . 
We use the fact that for a Gaussian random variable x 
with variance a 2 we have (e lx ) = e 1 ^^" / 2 . In this case, 
x = <jj m (t' — t) and <t 2 = 2D m (t — t'). We obtain after 
performing the integration 



N 



On — — / 

rn—1 



A h 



d„ 



lU> r , 



(C8) 
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